The following script and analyses were written by Colleen Bove and James Umbanhowar. This version of the script was last committed to GitHub on 22 Febuary 2020
In total, eleven studies met the standards of this meta-analysis, including the responses of thirteen Caribbean coral species collected from five different countries across the wider Caribbean (Figure 1 and Supplementary Table S1). Of the studies selected, only four performed fully factorial ocean acidification and warming experiments, and one performed two independent acidification and warming experiments. The most studied coral species from the Caribbean region were Porites astreoides (5 studies), Acropora cervicornis (4 studies), and Siderastrea siderea (4 studies). Finally, the Florida Keys and Belize were the two most-studied regions within the wider Caribbean fitting the criteria of this meta-analysis.
Meta-analysis of the dataset revealed that calcification rates of Caribbean corals were reduced by ocean warming but not ocean acidification (Figure 2 and Supplementary Figures S2). However, the 95% confidence interval of the combination of ocean warming and acidification overlapped zero, indicating no statistically clear trend towards synergistic or antagonistic effects of these treatments (Figure 2 and Supplementary Table S2).
Figure 1 Coral collection sites of all included studies with experimental study represented by colour and treatments represented by shape: acidification only (circle), warming only (triangle), and the combination of acidification and warming (square). The lower left insert displays close-up of the Belize collection sites and the upper right insert displays the Florida Keys collection sites.
# calculate ES and variance for ocean acidification
ES.oa <- escalc(measure = "SMDH", m1i = cm1, m2i = am3, sd1i = s1, sd2i = s3,
n1i = n1, n2i = n3, data = df)
# calculate ES and variance ES for ocean warming
ES.ow <- escalc(measure = "SMDH", m1i = cm1, m2i = wm4, sd1i = s1, sd2i = s4,
n1i = n1, n2i = n4, data = df)
ES.ow <- ES.ow %>% rename(yi.w = yi, vi.w = vi)
# create DF of calculated variables of OA
acid <- data.frame(cm1 = ES.oa$cm1, yi.a = ES.oa$yi, vi.a = ES.oa$vi)
# merge together warming and acid DF
ES.full <- merge(ES.ow, acid, by = "cm1")
ES.full <- ES.full[order(ES.full$cite), ]
ES.full$yi.w <- (ES.full$yi.w) * -1
ES.full$yi.a <- (ES.full$yi.a) * -1
### Standard deviation
# calculate sd for all
ES.full$sd1 <- ES.full$s1 * sqrt(ES.full$n1)
ES.full$sd3 <- ES.full$s3 * sqrt(ES.full$n3)
ES.full$sd4 <- ES.full$s4 * sqrt(ES.full$n4)
ES.full$sd6 <- ES.full$s6 * sqrt(ES.full$n6)
# calculate S for all treatments
ES.full$Sa <- sqrt((ES.full$sd1^2/ES.full$n1) + (ES.full$sd3^2/ES.full$n3))
ES.full$Sw <- sqrt((ES.full$sd1^2/ES.full$n1) + (ES.full$sd4^2/ES.full$n4))
ES.full$Si <- sqrt((ES.full$sd1^2/ES.full$n1) + (ES.full$sd6^2/ES.full$n6))
### Effect Size of Interaction term
# calculate ES using Harvey et al equations for acidification, warming, and
# interaction studies
ES.full$lnRRa <- ((ES.full$am3 - ES.full$cm1)/ES.full$Sa)
ES.full$lnRRt <- ((ES.full$wm4 - ES.full$cm1)/ES.full$Sw)
ES.full$yi.i <- ((ES.full$m6 - ES.full$wm4) - (ES.full$am3 - ES.full$cm1))/(2 *
ES.full$Si)
# calculate V for interaction studies ONLY
ES.full$vi.i <- (1/ES.full$n1) + (1/ES.full$n3) + (1/ES.full$n4) + (1/ES.full$n6) +
((ES.full$yi.i^2)/(2 * (ES.full$n1 + ES.full$n3 + ES.full$n4 + ES.full$n6)))
write.csv(ES.full, file = "~/Bove2020_MetaAnalysis_FrontMarSci/Files/ES.full.csv")
# read in the saved CSV and then removed any unneeded columns
df<-read.csv("~/Bove2020_MetaAnalysis_FrontMarSci/Files/ES.full.csv")
df<-df[,-c(2, 12:13, 16:18, 20:22, 24:26, 41:49)] #remove any unneeded columns
### data manipulation
#subset df for variance (v) values only as new df
df2b<-subset(df, select=c(X,vi.i,vi.a,vi.w))
#subset df to remove only variance (v) values as another new df
df.m <-subset(df, select=-c(vi.i,vi.a,vi.w))
# this creates df with treat and V as column headers in long format
df2b<-gather(df2b, treat, v,vi.i:vi.w)
#substitute y for v for merge later
df2b$treat<-gsub('v','y', df2b$treat)
#creates another long format df with treat and Y as headers (of Y values)
df.m <-gather(df.m, treat, y,yi.w:yi.i)
df.m $treat<-as.factor(df.m $treat)
#merge two df into one
df2<-merge(df.m, df2b, by=c("X", "treat"))
#remove any NAs
#df2<-na.omit(df2)
df2$study <- paste(df2$author, df2$year, sep = "_")
df2$X<-NULL
## centering treatments
df2$mag.p <- df2$a.pco2 - df2$c.pco2
df2$mag.t <- df2$w.temp - df2$c.temp
df2$mag.cent.p <- df2$mag.p - 127 # make a centered around 0 column
df2$mag.cent.t <- df2$mag.t - 1.0 # make a centered around 0 column
# subset data by treatment
OW<-subset(df2, treat=="yi.w")
OA<-subset(df2, treat=="yi.a")
IN<-subset(df2, treat=="yi.i")
# fit model to each treatment
ow <- rma(y, v, data=OW)
oa <- rma(y, v, data=OA)
both <- rma(y, v, data=IN)
# wide to long format df
df2 <- gather(data = df2, key = "treat.mag", value = "mag", mag.p:mag.t)
# remove NA values from df
completeFun <- function(data, desiredCols) {
completeVec <- complete.cases(data[, desiredCols])
return(data[completeVec, ])
}
df2 <- completeFun(df2, "y")
df2 <- completeFun(df2, "mag")
global.mod <- rma.mv(y,v, mods = ~ treat , random = list(~ (treat|study), ~ (treat|spec)), data = df2)
preds <- as.data.frame(predict(global.mod, newdata = df2))[,c(1,3:4)] # pull the mean and CI with the predict function
global.out <- cbind(df2, preds) # add columns to existing df for global output
#global.mod
global.out$treat <- factor(global.out$treat, levels = c("yi.a","yi.w","yi.i"))
global.out$wt <- 1/(global.out$v)
global_CI <-
ggplot()+
theme(panel.grid.major=element_blank(), legend.position="none", panel.grid.minor=element_blank(), panel.background=element_blank())+
theme(axis.line=element_line(color="black"))+
theme(legend.key=element_blank())+
ylab("Mean effect size (SMD)")+
xlab("")+
#ggtitle("Figure 2. Global standardized mean difference by treatment")+
scale_x_discrete(labels=c("yi.a" = "acidification", "yi.w" = "warming", "yi.i" = "combination"))+
geom_hline(yintercept=0, linetype=2, colour="black")+
geom_point(aes(x = treat, y = y, size = wt), colour="grey", data = global.out, alpha = 0.5, position=position_jitter(width=0.1))+
#geom_errorbar(aes(ymin=lb, ymax=ub, width=0.1), size=0.6)+
geom_linerange(data = global.out, aes(x = treat, ymin=ci.lb, ymax=ci.ub), size=0.8)+
#ylim(-1.6,0.6)+
geom_point(data = global.out, aes(x = treat, y = pred), size=10, shape=95)
div(ggplotly(global_CI), align = "center")
ggsave("~/Bove2020_MetaAnalysis_FrontMarSci/Figures/Manuscript/Figure2_GlobalCI.pdf", width = 5, height = 3.5)
Figure 2 Mean effect (standard mean difference) and 95% confidence interval of ocean acidification, warming, and the combination of acidification and warming on calcification rate for all studies in the meta-analysis. Grey circles indicate the effect size of each individual study and the size of each circle represents the weight of each study (1/SE). Clear statistical evidence of a treatment effect is identified when the 95% confidence interval does not overlap zero.
Corals from Belize only exhibited clearly reduced calcification rates under ocean warming (Figure 3A and Supplementary Table 3), while acidification, warming, and the combination of both stressors did not clearly alter experimental calcification rates of corals from the Florida Keys (Figure 3B and Supplementary Table S3). Further, the resulting QE suggests there is significant between-study variation (Supplementary Tables S3).
df2a<-subset(df2, region!="Curacao")
df2a<-subset(df2a, region!="Little Cayman")
df2a<-subset(df2a, region!="Mexico")
region.mod <- rma.mv(y,v, mods = ~ treat * region, random = list(~ (treat|study), ~ (treat|spec)), data = df2a)
preds2 <- as.data.frame(predict(region.mod, newdata = df2a))[,c(1,3:4)] # pull the mean and CI with the predict function
region.out <- cbind(df2a, preds2) # add columns to existing df for global output
#region.mod
region.out$treat <- factor(region.out$treat, levels = c("yi.a","yi.w","yi.i"))
region.out$wt <- 1/(region.out$v)
region_CI <- ggplot()+
theme(panel.grid.major=element_blank(), legend.position="none", panel.grid.minor=element_blank(), panel.background=element_blank())+
theme(axis.line=element_line(color="black"))+
theme(legend.key=element_blank())+
ylab("Mean effect size (SMD)")+
xlab("")+
scale_x_discrete(labels=c("yi.a" = "acidification", "yi.w" = "warming", "yi.i" = "combination"))+
geom_hline(yintercept=0, linetype=2, colour="black")+
geom_point(aes(x = treat, y = y, size = wt), colour="grey", data = region.out, alpha = 0.5, position=position_jitter(width=0.1))+
#ggtitle("Figure 3.Standardized mean difference by treatment for Florida Keys and Southern Belize")+
#geom_errorbar(aes(ymin=lb, ymax=ub, width=0.1), size=0.6)+
geom_linerange(data = region.out, aes(x = treat, ymin=ci.lb, ymax=ci.ub), size=0.8)+
geom_point(data = region.out, aes(x = treat, y = pred), size=10, shape=95)+
facet_wrap(~region)
div(ggplotly(region_CI), align = "center")
ggsave("~/Bove2020_MetaAnalysis_FrontMarSci/Figures/Manuscript/Figure3_RegionalCI.pdf", width = 8, height = 3.5)
Figure 3 Mean effect (standard mean difference) and 95% confidence interval of ocean acidification, warming, and the combination of acidification and warming on calcification rate for (A) Florida Keys corals and (B) Belize corals. Grey circles indicate the effect size of each individual study and the size of each circle represents the weight of each study (1/SE). Clear statistical evidence of a treatment effect is identified when the 95% confidence interval does not overlap zero.
Secondary analysis of mean calcification rates (mg cm−2 day−1) against treatment temperature across all Florida and Belize studies revealed a parabolic response to temperature (Figure 4 and Supplementary Tables S4, S5). Similarly, mean calcification rates across ΩArag resulted in a nonlinear response to acidification (Figure 5 and Supplementary Tables S6, S7). Both nonlinear trends in response to temperature and ΩArag were a result of treatment rather than region (Supplementary Tables S6, S7), suggesting regional differences identified in the meta-analysis were due to experimental designs employed to represent current regional environmental differences.
calc.w<-calc[,-c(15:18, 23:27)] #remove columns not needed for warming analysis
# remove and NA columns from warming column
calc.w <- completeFun(calc.w, "wm4")
# remove mexico studies
calc.w<-subset(calc.w, region != "Mexico")
# subset calc for calcication rates
calc.b<-subset(calc.w, select=c(num, cm1, wm4))
# create df with treat and rate as column headers in long format
calc.b<-gather(calc.b, treat, rate, cm1:wm4)
# rename for useful treatment names
calc.b$treat<-gsub('cm1','control', calc.b$treat)
calc.b$treat<-gsub('wm4','warm', calc.b$treat)
#subset calc for SE
calc.c<-subset(calc.w, select=c(num, s1, s4))
# create df with treat and SE as column headers in long format
calc.c<-gather(calc.c, treat, SE, s1:s4)
# rename for useful treatment names
calc.c$treat<-gsub('s1','control', calc.c$treat)
calc.c$treat<-gsub('s4','warm', calc.c$treat)
#subset calc to remove only columns from above
calc.d <-subset(calc.w, select=-c(s1 ,s4 ,cm1, wm4, c.pco2, n1, n4))
# this creates df with treat and V as column headers in long format
calc.d<-gather(calc.d, treat, temp, c.temp:w.temp)
# rename for useful treatment names
calc.d$treat<-gsub('c.temp','control', calc.d$treat)
calc.d$treat<-gsub('w.temp','warm', calc.d$treat)
#merge dfs into one, final warming
calc.m<-merge(calc.d, calc.b, by=c("num", "treat"))
calc.warm<-merge(calc.m, calc.c, by=c("num", "treat"))
# create column for weight
calc.warm$wt <- 1 / calc.warm$SE
# create study column
calc.warm$study <- paste(calc.warm$author, calc.warm$year, sep = "_")
# center temps
calc.warm$cent.temp <- calc.warm$temp - mean(calc.warm$temp)
calc.warm$cent.pco2 <- rep("NA", 60)
calc.warm$scale.temp <- scale(calc.warm$temp, center = TRUE, scale = TRUE)
warm.rate.region.mod <- lmer(rate ~ scale.temp + region + days + (1|study) + (1|spec), weights = 1/SE, data = calc.warm,REML=F) #sing fit
warm.rate.region2.mod <- lmer(rate ~ scale.temp + region + I(scale.temp^2) + days + (1|study) + (1|spec), weights = 1/SE, data = calc.warm,REML=F)
warm.rate2.mod <- lmer(rate ~ scale.temp + I(scale.temp^2) + days + (1|study) + (1|spec), weights = 1/SE, data = calc.warm,REML=F)
warm.rate.mod <- lmer(rate ~ scale.temp + days+ (1|study) + (1|spec), weights = 1/SE, data = calc.warm,REML=F)
warm.region.mod <- lmer(rate ~ region + days + (1|study) + (1|spec), weights = 1/SE, data = calc.warm,REML=F) #sing fit
#summary(warm.rate.region2.mod) # view summary output of model
#confint(warm.rate.region2.mod, method="boot")
warm.pred <- data.frame(temp=seq(26,32.01,length.out=100), region=rep(c("Belize", "Florida Keys"), length.out = 100), pred=predict(warm.rate.region2.mod, newdata = data.frame(scale.temp=seq(min(calc.warm$scale.temp),max(calc.warm$scale.temp),length.out=100), days=seq(min(calc.warm$days),max(calc.warm$days),length.out=100), region=rep(c("Belize", "Florida Keys"), length.out = 100)), re.form=NA))
global_OW_calc <-
ggplot(calc.warm, aes(x=temp, y=rate))+
geom_point(aes(size=1/SE, shape=region, fill=region), colour="black")+
theme(panel.grid.major=element_blank(), legend.position="right", panel.grid.minor=element_blank(), panel.background=element_blank())+
theme(axis.line=element_line(color="black"))+
theme(legend.key=element_blank())+
scale_shape_manual("Region", values = c(21,22,24))+
scale_fill_manual("Region", values = c("#5ab4ac", "#d8b365", "#f5f5f5"))+
scale_color_manual("Region", values = c("#5ab4ac", "#d8b365", "#f5f5f5"))+
geom_line(data = warm.pred, aes(y = pred, group = region, colour = region), size = 1)+
ylab(bquote('Calcification rate (mg' ~cm^-2* ~day^-1*')'))+
xlab(bquote('Temperature ('*degree*C*')'))
global_OW_calc
ggsave("~/Bove2020_MetaAnalysis_FrontMarSci/Figures/Manuscript/Figure4_CalcWarm.pdf", width = 7, height = 4)
Figure 4. Mean calcification rate (mg cm−2 day−1) of corals from each warming study by treatment temperature (°C) with the linear mixed effects model quadratic fit. Shape and colour of each point denotes study region (blue circle = Belize; brown square = Florida) and size of shape represents the weight (1/SE) of study.
calc.a<-calc[,-c(19:27)] #remove columns not needed for acid analysis
calc.a$p.mag <- calc.a$c.pco2 - calc.a$a.pco2 # pCO2 magnitude
calc.a <- completeFun(calc.a, "am3") # subset only for acidification studies
calc.a <- subset(calc.a, region != "Curacao")
#subset calc for calcication rates
calc.b<-subset(calc.a, select=c(num, cm1, am3))
# create df with treat and rate as column headers in long format
calc.b<-gather(calc.b, treat, rate, cm1:am3)
# rename for useful treatment names
calc.b$treat<-gsub('cm1','control', calc.b$treat)
calc.b$treat<-gsub('am3','acid', calc.b$treat)
#subset calc for SE
calc.c<-subset(calc.a, select=c(num, s1, s3))
# create df with treat and SE as column headers in long format
calc.c<-gather(calc.c, treat, SE, s1:s3)
# rename for useful treatment names
calc.c$treat<-gsub('s1','control', calc.c$treat)
calc.c$treat<-gsub('s3','acid', calc.c$treat)
#subset calc to remove only columns from above
calc.d <-subset(calc.a, select=-c(s1 ,s3 ,cm1, am3, c.temp, n1, n3))
# this creates df with treat and V as column headers in long format
calc.d<-gather(calc.d, treat, acid, c.pco2:a.pco2)
# rename for useful treatment names
calc.d$treat<-gsub('c.pco2','control', calc.d$treat)
calc.d$treat<-gsub('a.pco2','acid', calc.d$treat)
#merge dfs into one, final acid
calc.m<-merge(calc.d, calc.b, by=c("num", "treat"))
calc.arag<-merge(calc.m, calc.c, by=c("num", "treat"))
# create column for weight
calc.arag$wt <- 1 / calc.arag$SE
# create study column
calc.arag$study <- paste(calc.arag$author, calc.arag$year, sep = "_")
# center temps
calc.arag$cent.pco2 <- calc.arag$acid - mean(calc.arag$acid)
#calc.arag$cent.pco2 <- rep("NA", 66)
calc.arag$scale.arag <- scale(calc.arag$acid, center = TRUE, scale = TRUE)
acid.rate.region.mod <- lmer(rate ~ scale.arag + region + days + (1|study) + (1|spec), weights = 1/SE, data = calc.arag,REML=F) #sing fit
acid.rate.region2.mod <- lmer(rate ~ scale.arag + region + I(scale.arag^2) + days + (1|study) + (1|spec), weights = 1/SE, data = calc.arag,REML=F)
acid.rate2.mod <- lmer(rate ~ scale.arag + I(scale.arag^2) + days + (1|study) + (1|spec), weights = 1/SE, data = calc.arag,REML=F)
acid.rate.mod <- lmer(rate ~ scale.arag + days+ (1|study) + (1|spec), weights = 1/SE, data = calc.arag,REML=F)
acid.region.mod <- lmer(rate ~ region + days + (1|study) + (1|spec), weights = 1/SE, data = calc.arag,REML=F) #sing fit
#summary(acid.rate.region2.mod) # continuing with this one for simplicity
#confint(acid.rate.region2.mod, method="boot")
acid.pred <- data.frame(acid=seq(0.7,5.8,length.out=100), region=rep(c("Belize", "Florida Keys"), length.out = 100),pred=predict(acid.rate.region2.mod, newdata = data.frame(scale.arag=seq(min(calc.arag$scale.arag),max(calc.arag$scale.arag),length.out=100), days=seq(min(calc.arag$days),max(calc.arag$days), length.out=100), region=rep(c("Belize", "Florida Keys"), length.out=100)),re.form=NA))
global_OA_calc <-
ggplot(calc.arag, aes(x=acid, y=rate))+
geom_point(aes(size=1/SE, shape=region, fill=region), colour="black")+
theme(panel.grid.major=element_blank(), legend.position="right", panel.grid.minor=element_blank(), panel.background=element_blank())+
theme(axis.line=element_line(color="black"))+
theme(legend.key=element_blank())+
scale_shape_manual("Region", values = c(21,22,24))+
scale_fill_manual("Region", values = c("#5ab4ac", "#d8b365", "#f5f5f5"))+
scale_color_manual("Region", values = c("#5ab4ac", "#d8b365", "#f5f5f5"))+
geom_line(data = acid.pred, aes(y = pred, group = region, colour = region), size = 1)+
ylab(bquote('Calcification rate (mg' ~cm^-2* ~day^-1*')'))+
xlab(bquote('Aragonite Saturation State'))
global_OA_calc
ggsave("~/Bove2020_MetaAnalysis_FrontMarSci/Figures/Manuscript/Figure5_CalcAcid.pdf", width = 7, height = 4)
Figure 5. Mean calcification rate (mg cm−2 day−1) of corals from each acidification study by treatment aragonite saturation state (ΩArag) with the linear mixed effects model quadratic fit. Shape and colour of each point denotes study region (blue circle = Belize; brown square = Florida Keys) and size of shape represents the weight (1/SE) of study.
Quantification of experimental design parameters within warming studies identified that magnitude of treatment, irradiance, seawater type used, and feeding frequency all clearly impacted calcification rates (Table 1). Specifically, studies that utilized natural seawater and those with a larger difference between the control and treatment temperatures within a study exhibited higher effect sizes, suggesting a less negative effect of treatment. Studies that employed higher irradiance levels in their systems demonstrated more negative effects of treatment light level. Finally, studies that reported feeding their corals twice a week were less impacted by warming treatment than those feeding three times a week, however, studies with no data on feeding were the least affected by treatment. Duration of experiment was deemed redundant in the model and was thus dropped.
Within the acidification studies, irradiance, seawater type used, feeding frequency, duration, and the interaction of duration with treatment magnitude impacted effect sizes, while magnitude alone was not clearly different (Table 1). Studies using natural seawater, employing higher irradiance levels, and those with longer durations resulted in great effect sizes, suggesting they lessened the effects of acidification treatment on calcification responses. Similar to warming studies, acidification studies in which feeding corals was conducted twice a week exhibited less negative responses to treatment than those feeding three times a week, with studies reporting no feeding data exhibiting the least negative responses to treatment. Finally, coral calcification responses were less impacted by acidification in studies with longer duration of exposure and a greater pCO2 change.
Table 1. Effect size estimate and 95% confidence intervals of experimental design parameters calculated for all wamring and acidification experiments.
full.param <- rbind(w.global.out[,-4], a.global.out[,-4])
param.tab <- kable(full.param, digits = 3) %>%
kable_styling(font_size = 14, full_width = F) %>%
pack_rows("Warming Experiments", 1, 6) %>%
pack_rows("Acidification Experiments", 7, 13)
param.tab
| Estimate | Lower 95% CI | Upper 95% CI | col | |
|---|---|---|---|---|
| Warming Experiments | ||||
| Intercept | -29.392 | -39.564 | -19.221 | negative |
| Magnitude | 10.341 | 5.580 | 15.101 | positive |
| Irradiance | -0.237 | -0.336 | -0.137 | negative |
| Seawater | 29.767 | 19.719 | 39.815 | positive |
| Feeding (3x) | -3.544 | -4.967 | -2.120 | negative |
| Feeding (N.D.) | 87.160 | 49.685 | 124.635 | positive |
| Acidification Experiments | ||||
| Intercept | -2.033 | -2.785 | -1.281 | negative |
| Magnitude | 0.000 | -0.002 | 0.003 | positive |
| Irradiance | 0.009 | 0.004 | 0.014 | positive |
| Seawater | 2.392 | 1.144 | 3.640 | positive |
| Feeding (3x) | -15.042 | -22.607 | -7.477 | negative |
| Feeding (N.D.) | 3.611 | 0.280 | 6.942 | positive |
| Duration | 0.206 | 0.090 | 0.323 | positive |
| Magnitude x Duration | -0.001 | -0.001 | 0.000 | negative |
Supplementary Table S1. List of included studies in meta-analysis with experimental design parameters. Treatment temperature and pCO2 reported are for treatment levels used in this analysis. Missing data from studies is denoted by N.D.
table1 <- kable(meta_tab2, booktabs = T) %>%
kable_styling(font_size = 10, full_width = F) %>%
column_spec(5, italic = T)
table1
| Study | Location | Collection.Sites | Species | Treatments | Experiment.Duration | Feeding | Seawater | Control.Temperature | Temperature | Control.pCO2 | pCO2 | irradiance | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | Albright et al. 2011 | Florida | Little Grecian, Key Largo | Porites astreoides | pCO2 | 49 d | NA | Natural seawater | 28.0 | NA | 380 | 800.0 | NA |
| 2 | Bedwell-Ivers et al. 2017 | Little Cayman | Little Cayman Research Center | Acropora cervicornis Porites divaricata | pCO2 | 28 d | NA | Natural seawater | 30.0 | NA | 487 | 791.0 | 800 |
| 3 | Bove et al. 2019 | Belize | Inshore, Offshore (Sapadilla Cayes) | Porites astreoides, Pseudodiploria strigosa, Siderastrea siderea, Undaria tenuifolia | Warming, pCO2, combination | 93 d | Frozen and newly hatched Artemia every other day | Filtered non-native seawater | 28.0 | 31.0 | 405 | 701.0 | 300 |
| 4 | Castillo et al. 2014 | Belize | Nearshore, Backreef, Forereef (Sapadilla Cayes) | Siderastrea siderea | Warming, pCO2 | 95 d | Frozen Artemia every other day | Artificial seawater | 28.0 | 32.0 | 472 | 704.0 | 250 |
| 5 | Enochs et al. 2014 | Florida | Dade County | Acropora cervicornis | pCO2 | 42 d | NA | Natural seawater | 28.0 | NA | 513 | 870.9 | 178 |
| 6 | Grottoli et al. 2014 | Mexico | Puerto Morelos | Porites divaricata, Porites astreoides, Orbicella faveolata | Warming | 15 d | Fed | Natural seawater | 30.6 | 31.6 | NA | NA | 600 |
| 7 | Horvath et al. 2016 | Belize | Nearshore, Backreef, Forereef (Sapadilla Cayes) | Siderastrea siderea | Warming, pCO2, combination | 60 d | Frozen Artemia twice weekly | Artificial seawater | 28.0 | 31.0 | 426 | 890.0 | 250 |
| 8 | Jury et al. 2010 | Curacao | Seaquarium | Madracis auretenra | pCO2 | 2 h | Artemia nauplii twice weekly | Natural seawater | 28.0 | NA | 391 | 1480.0 | 200 |
| 9 | Okazaki et al. 2017 | Florida | Key Biscayne, Biscayne National Park, Key West | Acropora cervicornis, Agaricia agaricites, Dichocoenia stokesii, Montastraea cavernosa, Orbicella faveolata, Porites astreoides, Porites divaricata, Pseudodiploria strigosa, Siderastrea siderea, Siderastrea radians, Solenastrea hyades | Warming, pCO2, combination | 6 w | Live rotifers and larval twice weekly | Natural seawater | 27.0 | 30.3 | 400 | 1340.0 | 327 |
| 11 | Towle et al. 2015 | Florida | South Florida | Acropora cervicornis | Warming, pCO2, combination | 8 w | Dried zooplankton twice weekly | Natural seawater | 26.0 | 30.0 | 390 | 800.0 | 353 |
| 12 | Kenkel et al. 2013 | Florida | Inshore and Offshore Sugarloaf Key | Porites astreoides | Warming | 43 d | NA | Natural seawater | 27.2 | 30.9 | NA | NA | NA |
#save_kable(table1, "Table1_Studies_16Jan20.pdf")
Supplementary Table S2. Global meta-analysis mixed effects model (function rma.mv) output by treatment only, with random effects of study and species, used in Figure 2. The test for residual heterogeneity and significance is represented by QE (P-value) and the test of moderators is QM (P-value).
tidy.rma <- function(x, ...) {
with(
x,
data.frame(
estimate = b,
std.error = se,
z.value = zval,
p.value = pval,
conf.low = ci.lb,
conf.high = ci.ub
)
)
}
tidy.rma.Q <- function(x, ...) {
with(
x,
data.frame(
QE = QE,
p.value = QEp,
QM = QM,
p.value = QMp
)
)
}
global.Q <- tidy.rma.Q(global.mod)
global.Q <- data.frame(t(global.Q))
row.names(global.Q) <- c("QE", "p value", "QM", " p value")
colnames(global.Q) <- "estimate"
global.Q$std.error <- rep("", 4)
global.Q$z.value <- rep("", 4)
global.Q$p.value <- rep("", 4)
global.Q$conf.low <- rep("", 4)
global.Q$conf.high <- rep("", 4)
global.mod.tab <- tidy.rma(global.mod)
#row.names(global.mod.tab) <- c("acidification", "combination", "warming")
global.mod.tab$std.error <- round(global.mod.tab$std.error, 3)
global.mod.tab$z.value <- round(global.mod.tab$z.value, 3)
global.mod.tab$p.value <- round(global.mod.tab$p.value, 3)
global.mod.tab$conf.low <- round(global.mod.tab$conf.low, 3)
global.mod.tab$conf.high <- round(global.mod.tab$conf.high, 3)
global.mod.tab <- rbind(global.mod.tab, global.Q)
global.mod.tab <- global.mod.tab[,1:4]
table2 <- kable(global.mod.tab, digits = 2) %>%
kable_styling(font_size = 14, full_width = F) %>%
pack_rows("Model Results", 1, 3) %>%
pack_rows("Variance Components", 4, 7)
table2
| estimate | std.error | z.value | p.value | |
|---|---|---|---|---|
| Model Results | ||||
| intrcpt | -1.20 | 0.776 | -1.545 | 0.122 |
| treatyi.i | 1.19 | 1.41 | 0.843 | 0.399 |
| treatyi.w | -1.07 | 1.19 | -0.898 | 0.369 |
| Variance Components | ||||
| QE | 823.95 | |||
| p value | 0.00 | |||
| QM | 2.40 | |||
| p value | 0.30 | |||
#save_kable(table2, "Table2_GlobalMod_16Jan20.pdf")
Supplementary Table S3. Meta-analysis mixed effects model (function rma.mv) output of treatment by region (Belize versus Florida Keys), with random effects of study and species, used in Figure 3. The test for residual heterogeneity and significance is represented by QE (P-value) and the test of moderators is QM (P-value).
regional.Q <- tidy.rma.Q(region.mod)
regional.Q <- data.frame(t(regional.Q))
row.names(regional.Q) <- c("QE", "p value", "QM", " p value")
colnames(regional.Q) <- "estimate"
regional.Q$std.error <- rep("", 4)
regional.Q$z.value <- rep("", 4)
regional.Q$p.value <- rep("", 4)
regional.Q$conf.low <- rep("", 4)
regional.Q$conf.high <- rep("", 4)
region.mod.tab <- tidy.rma(region.mod)
#row.names(region.mod.tab) <- c("acidification", "combination", "warming","acidification ", "combination ", "warming ")
region.mod.tab$std.error <- round(region.mod.tab$std.error, 3)
region.mod.tab$z.value <- round(region.mod.tab$z.value, 3)
region.mod.tab$p.value <- round(region.mod.tab$p.value, 3)
region.mod.tab$conf.low <- round(region.mod.tab$conf.low, 3)
region.mod.tab$conf.high <- round(region.mod.tab$conf.high, 3)
region.mod.tab <- rbind(region.mod.tab, regional.Q)
region.mod.tab <- region.mod.tab[,1:4]
table3 <- kable(region.mod.tab, digits = 4, booktabs = T) %>%
kable_styling(font_size = 14, full_width = F) %>%
pack_rows("Belize", 1, 3) %>%
pack_rows("Florida", 4, 6) %>%
pack_rows("Variance Components", 7, 10)
table3
| estimate | std.error | z.value | p.value | |
|---|---|---|---|---|
| Belize | ||||
| intrcpt | -0.9288 | 1.324 | -0.701 | 0.483 |
| treatyi.i | 1.5625 | 2.24 | 0.697 | 0.486 |
| treatyi.w | -3.6439 | 2.069 | -1.761 | 0.078 |
| Florida | ||||
| regionFlorida Keys | 0.0861 | 1.879 | 0.046 | 0.963 |
| treatyi.i:regionFlorida Keys | -0.9606 | 3.184 | -0.302 | 0.763 |
| treatyi.w:regionFlorida Keys | 5.1661 | 3.061 | 1.688 | 0.091 |
| Variance Components | ||||
| QE | 504.3862 | |||
| p value | 0.0000 | |||
| QM | 9.3736 | |||
| p value | 0.0951 | |||
#save_kable(table3, "Table3_RegionalMod_16Jan20.pdf")
Supplementary Table S4. Temperature and region linear mixed effects model section using AICc. All models were run with random effects for study and species and a weight of 1/SE. The asterisk (*) denotes the selected model run for the final analysis.
warm.mag.tab <- AICc(warm.rate.region.mod, warm.rate.region2.mod, warm.rate2.mod, warm.rate.mod, warm.region.mod)
tableS4 <- kable(warm.mag.tab, digits = 2, booktabs = T) %>%
kable_styling(font_size = 14, full_width = F) %>%
row_spec(2, bold = T)
tableS4
| df | AICc | |
|---|---|---|
| warm.rate.region.mod | 7 | 109.99 |
| warm.rate.region2.mod | 8 | 75.39 |
| warm.rate2.mod | 7 | 72.85 |
| warm.rate.mod | 6 | 107.84 |
| warm.region.mod | 6 | 125.18 |
Supplementary Table S5. Temperature and region best fit linear mixed effects model output and 95% confidence intervals of the calcification rates in response to treatment temperature plotted in Figure 4.
tab_model(warm.rate.region2.mod, show.r2 = FALSE, show.aicc = TRUE)
| Â | rate | ||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 0.74 | -1.72 – 3.21 | 0.555 |
| scale.temp | -0.23 | -0.31 – -0.16 | <0.001 |
| region: Florida Keys | 0.24 | -1.05 – 1.53 | 0.713 |
| I(scale.temp^2) | -0.36 | -0.44 – -0.27 | <0.001 |
| days | 0.01 | -0.02 – 0.04 | 0.659 |
| Random Effects | |||
| σ2 | 0.28 | ||
| τ00 spec | 0.20 | ||
| τ00 study | 0.18 | ||
| ICC | 0.57 | ||
| N study | 5 | ||
| N spec | 11 | ||
| Observations | 60 | ||
| AIC | 75.386 | ||
Supplementary Table S6. Aragonite saturation state and region linear mixed effects model section using AICc. All models were run with random effects for study and species and a weight of 1/SE. The asterisk (*) denotes the bet fit model.
arag.mag.tab <- AICc(acid.rate.region.mod, acid.rate.region2.mod, acid.rate2.mod, acid.rate.mod, acid.region.mod)
tableS6 <- kable(arag.mag.tab, digits = 2, booktabs = T) %>%
kable_styling(font_size = 14, full_width = F) %>%
row_spec(2, bold = T)
tableS6
| df | AICc | |
|---|---|---|
| acid.rate.region.mod | 7 | 57.97 |
| acid.rate.region2.mod | 8 | 49.85 |
| acid.rate2.mod | 7 | 47.98 |
| acid.rate.mod | 6 | 56.39 |
| acid.region.mod | 6 | 105.28 |
Supplementary Table S7. Aragonite saturation state (ΩArag) and region best fit linear mixed effects model output and 95% confidence intervals of the calcification rates in response to treatment ΩArag plotted in Figure 5.
tab_model(acid.rate.region2.mod, show.r2 = FALSE, show.aicc = TRUE)
| Â | rate | ||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | -0.19 | -3.06 – 2.69 | 0.900 |
| scale.arag | 0.10 | 0.07 – 0.14 | <0.001 |
| region: Florida Keys | 0.47 | -1.00 – 1.95 | 0.532 |
| I(scale.arag^2) | -0.04 | -0.07 – -0.02 | 0.001 |
| days | 0.02 | -0.02 – 0.05 | 0.376 |
| Random Effects | |||
| σ2 | 0.20 | ||
| τ00 spec | 0.17 | ||
| τ00 study | 0.25 | ||
| ICC | 0.68 | ||
| N study | 5 | ||
| N spec | 12 | ||
| Observations | 140 | ||
| AIC | 49.853 | ||
oa.for <- forest(oa, slab=paste(OA$cite, OA$spec, OA$site, sep=" - "), order=order(OA$cite))
Supplementary Figure S1. Forest plot depicting individual standard mean difference (SMD) and confidence interval of each species-site-study combination included in the meta-analysis of ocean acidification only studies. Studies with confidence intervals not overlapping zero (dotted line) denote either a significant increase (positive values) or decrease (negative values) in calcification response under experimental acidification compared to the corresponding control treatment.
oa.fun <- funnel(oa, level=c(90, 95, 99), shade=c("white", "gray55", "gray75"), legend=TRUE)
Supplementary Figure S2. Funnel plot of the standard mean difference (SMD) against standard error of each species-site-study combination included in the meta-analysis of ocean acidification only studies. Dotted lines represent confidence intervals of all studies, background plot colour represents statistical significance of individual study, the black vertical line denotes the overall effect of all included ocean acidification studies.
pdf("~/Bove2020_MetaAnalysis_FrontMarSci/Figures/Supplemental/FigureS1_OAforest.pdf", width = 11, height = 12)
forest(oa, slab=paste(OA$cite, OA$spec, OA$site, sep=" - "), order=order(OA$cite))
dev.off()
pdf("~/Bove2020_MetaAnalysis_FrontMarSci/Figures/Supplemental/FigureS2_OAfunnel.pdf", width = 7, height = 6)
funnel(oa, level=c(90, 95, 99), shade=c("white", "gray55", "gray75"), legend=TRUE)
dev.off()
ow.for <- forest(ow, slab=paste(OW$cite, OW$spec, OW$site, sep=" - "), order=order(OW$cite))
Supplementary Figure S3. Forest plot depicting individual standard mean difference (SMD) and confidence interval of each species-site-study combination included in the meta-analysis of ocean warming only studies. Studies with confidence intervals not overlapping zero (dotted line) denote either a significant increase (positive values) or decrease (negative values) in calcification response under experimental warming compared to the corresponding control treatment.
ow.fun <- funnel(ow, level=c(90, 95, 99), shade=c("white", "gray55", "gray75"), legend=TRUE)
Supplementary Figure S4. Funnel plot of the standard mean difference (SMD) against standard error of each species-site-study combination included in the meta-analysis of ocean warming only studies. Dotted lines represent confidence intervals of all studies, background plot colour represents statistical significance of individual study, the black vertical line denotes the overall effect of all included ocean acidification studies.
pdf("~/Bove2020_MetaAnalysis_FrontMarSci/Figures/Supplemental/FigureS3_OWforest.pdf", width = 9, height = 10)
forest(ow, slab=paste(OW$cite, OW$spec, OW$site, sep=" - "), order=order(OW$cite))
dev.off()
pdf("~/Bove2020_MetaAnalysis_FrontMarSci/Figures/Supplemental/FigureS4_OWfunnel.pdf", width = 7, height = 5)
funnel(ow, level=c(90, 95, 99), shade=c("white", "gray55", "gray75"), legend="topleft")
dev.off()
com.for <- forest(both, slab=paste(IN$cite, IN$spec, IN$site, sep=" - "), order=order(IN$cite))
Supplementary Figure S5. Forest plot depicting individual standard mean difference (SMD) and confidence interval of each species-site-study combination included in the meta-analysis of the combination of ocean acidification and warming studies.
com.fun <- funnel(both, level=c(90, 95, 99), shade=c("white", "gray55", "gray75"), legend=TRUE)
Supplementary Figure S6. Funnel plot of the standard mean difference (SMD) against standard error of each species-site-study combination included in the meta-analysis of combined ocean acidification and warming only. Dotted lines represent confidence intervals of all studies, background plot colour represents statistical significance of individual study, the black vertical line denotes the overall effect of all included ocean acidification studies.
pdf("~/Bove2020_MetaAnalysis_FrontMarSci/Figures/Supplemental/FigureS5_bothforest.pdf", width = 9, height = 8)
forest(both, slab=paste(IN$cite, IN$spec, IN$site, sep=" - "), order=order(IN$cite))
dev.off()
pdf("~/Bove2020_MetaAnalysis_FrontMarSci/Figures/Supplemental/FigureS6_bothfunnel.pdf", width = 7, height = 6)
funnel(both, level=c(90, 95, 99), shade=c("white", "gray55", "gray75"), legend=TRUE)
dev.off()
sessionInfo()
## R version 3.5.2 (2018-12-20)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] cowplot_1.0.0 kableExtra_1.1.0 sjlabelled_1.1.1
## [4] sjmisc_2.8.2 sjPlot_2.7.2 plotly_4.9.0
## [7] shiny_1.4.0 Rmisc_1.5 lattice_0.20-38
## [10] MuMIn_1.43.6 RColorBrewer_1.1-2 tidystats_0.4
## [13] forcats_0.4.0 stringr_1.4.0 dplyr_0.8.3
## [16] purrr_0.3.3 readr_1.3.1 tibble_2.1.3
## [19] tidyverse_1.2.1 lme4_1.1-21 plyr_1.8.4
## [22] tidyr_1.0.0 ggplot2_3.2.1 MAd_0.8-2.1
## [25] metafor_2.1-0 Matrix_1.2-17
##
## loaded via a namespace (and not attached):
## [1] TH.data_1.0-10 minqa_1.2.4 colorspace_1.4-1
## [4] ellipsis_0.3.0 estimability_1.3 parameters_0.2.0
## [7] rstudioapi_0.10 ggrepel_0.8.1 mvtnorm_1.0-11
## [10] lubridate_1.7.4 xml2_1.2.2 codetools_0.2-16
## [13] splines_3.5.2 mnormt_1.5-5 knitr_1.25
## [16] zeallot_0.1.0 jsonlite_1.6 nloptr_1.2.1
## [19] ggeffects_0.12.0 Cairo_1.5-10 broom_0.5.2
## [22] compiler_3.5.2 httr_1.4.1 sjstats_0.17.6
## [25] emmeans_1.4.2 backports_1.1.5 assertthat_0.2.1
## [28] fastmap_1.0.1 lazyeval_0.2.2 cli_1.1.0
## [31] later_1.0.0 formatR_1.7 htmltools_0.4.0
## [34] tools_3.5.2 coda_0.19-3 gtable_0.3.0
## [37] glue_1.3.1 Rcpp_1.0.2 cellranger_1.1.0
## [40] vctrs_0.2.0 nlme_3.1-141 crosstalk_1.0.0
## [43] psych_1.8.12 insight_0.6.0 xfun_0.10
## [46] rvest_0.3.4 mime_0.7 lifecycle_0.1.0
## [49] MASS_7.3-51.4 zoo_1.8-6 scales_1.0.0
## [52] hms_0.5.1 promises_1.1.0 parallel_3.5.2
## [55] sandwich_2.5-1 yaml_2.2.0 stringi_1.4.3
## [58] highr_0.8 bayestestR_0.4.0 boot_1.3-23
## [61] rlang_0.4.1 pkgconfig_2.0.3 evaluate_0.14
## [64] htmlwidgets_1.5.1 labeling_0.3 tidyselect_0.2.5
## [67] magrittr_1.5 R6_2.4.0 generics_0.0.2
## [70] multcomp_1.4-10 pillar_1.4.2 haven_2.1.1
## [73] foreign_0.8-72 withr_2.1.2 survival_2.44-1.1
## [76] performance_0.4.0 modelr_0.1.5 crayon_1.3.4
## [79] rmarkdown_1.16 grid_3.5.2 readxl_1.3.1
## [82] data.table_1.12.6 digest_0.6.22 webshot_0.5.1
## [85] xtable_1.8-4 httpuv_1.5.2 stats4_3.5.2
## [88] munsell_0.5.0 viridisLite_0.3.0